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1  Introduction 


Modeling  of  small  thruster  expansions  has  seen  a  sustained  development  of  computational 
methods  and  options  over  the  last  decade.  The  physics  of  such  flows  encompasses  multiple 
length-scales  from  continuum  to  rarefied  as  the  gas  expands  through  the  nozzle  into  the 
near-space  vacuum.  Moreover,  with  the  advent  of  MEMS  fabrication  both  cold  and  hot 
gas  working  fluids  experience  non-continuum  effects  because  the  surface-to- volume  ratio  is 
higher  than  in  traditional  axi- symmetric  nozzles.  As  MEMS  fabrication  has  improved  to 
enable  higher  stagnation  pressures  particle  methods  are  still  relevant.  However,  they  must 
be  able  to  take  advantage  of  the  large  collision  rates  in  the  early  part  of  the  expansions, 
rather,  than  be  penalized  by  enormous  computational  costs  of  DSMC.  Finally,  increasing 
the  stagnation  temperature  does  improve  thrust  levels  (although  not  efficiency),  but, 
predicting  the  sustainability  of  the  nozzle  material  means  that  particle-flow  and  thermal 
simulations  must  be  coupled. 

Because  the  residence  time  of  the  gas  from  the  nozzle  throat  to  the  exit  is  sufficiently 
short  the  chemistry  of  multi-species  flows  through  meso  and  MEMS  -scaled  nozzles  is  es¬ 
sentially  frozen.  However,  given  the  rapidly  falling  gas  temperature  caused  by  the  strong 
supersonic  expansion  condensation  both  inside  and  beyond  the  nozzle  exit  are  likely  in 
many  neutral  gas  flows.  Whether  the  condensation  is  homogeneous  or  heterogeneous  is  an 
important  question.  From  an  operational  point  of  view,  heterogeneous  condensation  rates 
are  orders  of  magnitude  faster  than  homogeneous  condensation,  however,  scientifically 
such  processes  are  more  difficult  to  model  due  to  the  lack  of  detailed  information  re¬ 
garding  trace  impurities  which  initiate  the  process.  Condensation  is  an  important  source 
of  spacecraft  contamination  and  in  such  a  context  can  only  be  studied  using  particle 
approaches. 

The  research  that  has  been  conducted  in  my  group  has  sought  to  address  the  above 
technical  challenges  and  issues  and  will  be  summarized  in  these  lecture  notes.  It  should 
be  stated  from  the  out  set  that  the  modeling  of  MEMS  and  small  nozzle  thruster  simu¬ 
lations  using  particle  and  continuum  approaches  has  been  carried  out  by  others,  notably 
Ivanov, (1)  Gimelshein  and  Wysong,(2)  Aleexenko,(3)  and  others  and  the  experimental 
endeavors  to  provide  modelers  with  important  thrust  measurements  include  the  work  of 
Kestdever,(4;  5)  Hitt,  (6),  for  example.  Due  to  the  lack  of  time,  these  notable  contributions 
will  not  be  summarized  although  an  attempt  will  be  made  to  provide  appropriate  attribu¬ 
tion.  The  author  apologizes  in  advance  for  any  omissions.  The  organization  of  the  lecture 
is  as  follows.  After  discussing  the  numerous  applications  of  modeling  supersonic  expan¬ 
sions  for  finite  Knudsen  number  flows,  the  inter-relationships  among  the  different  types 
of  particle  computational  tools  used  in  the  research  here  is  discussed.  Then  the  physics 
of  two  important  types  of  expanding  flows  are  examined.  First  we  discuss  MEMS  flows 
for  simple  gases.  We  then  discuss  what  we  have  learned  from  our  modeling  of  condensing 
flows  covering  the  spectrum  of  chemically  inert  or  argon  flows  through  small  polyatomic 
systems  such  as  water,  carbon  dioxide,  ammonia,  and  ethanol.  Conclusions  in  terms  of 
anticipated  future  directions  and  challenges  are  then  presented. 
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2  Particle  methods  Computational  Tools 

2.1  DSMC  and  Collision- limiter  Techniques 

The  direct  simulation  Monte  Carlo  method  first  proposed  by  Bird(7)  is  frequently  de¬ 
scribed  as  a  phenomenological  approach  for  simulation  of  rarefied  gas  flows  by  modeling 
the  motion  of  fictitious  particles.  The  correct  usage  of  the  DSMC  method  requires  time 
and  space  discretization  in  sufficiently  small  quantities  to  ensure  that  the  predicted  flow 
parameters  correspond  to  a  solution  of  the  Boltzmann  equation.  The  latter  is  the  most 
general  equation  to  describe  the  transport  of  gases  in  a  fully  kinetic  manner.  DSMC  is  a 
numerical  approach  for  solving  the  Boltzmann  equation,  an  integro-differential  equation 
that  in  general  is  not  amenable  to  analytic  solution,  if  the  gas  is  dilute,  composed  of  binary 
collisions  and  molecular  chaos  holds.  As  has  been  discussed  by  Ivanov  and  Rogasinky(8) 
and  later  Alexeenko,(9)  the  phenomenological  description  of  the  method  is  not  neces¬ 
sary.  The  two  widely  used  approximate  numerical  schemes  for  the  collisional  relaxation 
used  in  the  direct  simulation  Monte  Carlo  method  are  the  No  Time  Counter  (NTC)  due 
to  Bird(7)  and  the  Majorant  Frequency  Scheme  (MFS)  of  Ivanov  and  Rogasinky.  Both 
schemes  have  the  same  numerical  efficiency,  with  computational  cost  linearly  proportional 
to  the  number  of  particles,  N  instead  of  N2.  The  MFS,  however,  is  more  accurate  than 
NTC  because  it  allows  one  to  model  the  collisional  relaxation  with  a  smaller  number  of 
particles. 

The  MFS  relaxation  scheme  during  a  time  step,  At  can  be  summarized  as  follows. 
The  initial  state  Co  is  selected  from  the  probability  density  f^(C,  0)  where  the  velocities 
of  N  particles  can  be  described  by, 

C  =  (vi,v2,  ....vN) 

and  /w  (Ci,  t)  is  the  the  probability  density  for  the  stochastic  vector  C  to  have  value  the 
value  C\  at  time  t.  Then  the  time  between  two  collisions,  rm ,  is  determined  by  setting 

rm  =  — z'”1  In  rand,  (1) 


where  vm  is  the  majorant  frequency.  The  time  counter  within  the  time  step  is  increased 
by 


tn+l  h  f  7~n 


and  if  tn+\  >  At  then  no  more  particles  are  selected  for  collisions  in  the  present  time 
step.  Else,  a  new  pair  of  particles  with  velocities  ( Vi,Vj )  are  selected  from  the  system  of 
N  particles  if, 


rand  < 


giMdij) 

\.9ij(J(9ij)\max 


where  g.,j  is  the  relative  collision  velocity  and  a  is  the  collision  cross  section.  If  the 
inequality  holds  then  the  selected  pair  is  allowed  to  collide  and  they  obtain  new  velocities 
after  the  collision  of  ( vd,Vjf ).  Else  one  determines  the  time  of  the  next  collision. 

It  should  be  noted  that  unlike  the  Navier-Stokes  equations,  the  Boltzmann  equation 
of  transport  is  valid  over  all  flow  regimes.  In  the  high-collision  limit,  or  the  very  low 
Knudsen  number  limit,  the  Boltzmann  equation  can  be  shown  using  Chapman-Enskog 
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theory  to  reduce  to  the  Navier-Stokes  and  Euler’s  equations.  Since  the  Knudsen  num¬ 
ber  may  be  related  to  the  Mach  and  the  Reynolds  numbers  the  selection  of  whether  the 
Navier-Stokes  or  the  Euler’s  inviscid  flow  equation  is  appropriate  depends  on  the  degree  of 
viscosity  in  the  flow.  In  particular,  if  the  flow  is  inviscid  the  Knudsen  number  approaches 
zero,  all  collisions  are  perfectly  elastic  and  reversible,  and  the  distribution  function  is  the 
equilibrium  Maxwell-Boltzmann  velocity  distribution  function.  While  this  is  an  idealiza¬ 
tion,  the  core  region  in  a  supersonic  nozzle  expansion  is  a  good  example  of  where  such  a 
situation  physically  occurs.  A  great  deal  of  computation  time  can  be  saved  if  one  knows  a 
priori  that  equilibrium  occurs  after  some  specified  number  of  collisions  (i.e.,  the  collision 
limiter  or  enforcer)  and  the  MFS  algorithm  can  be  naturally  adapted  to  implement  this 
procedure  within  a  time  step  within  a  DSMC  cell  as  follows.  Using  the  same  initial  state 
as  above,  select  particles  subject  to  the  inequality  given  by  Eq.  I.  If  the  selected  pair 
satisfies  the  inequality  perform  a  collision  and  assign  new  velocities  as  in  the  MFS  proce¬ 
dure.  If  not,  select  another  pair  and  use  Eq.  1  to  test  again  until  the  necessary  number 
of  collisions  (as  specified  by  the  collision  limiter  number)  have  been  performed.  When  so 
accomplished,  move  on  to  the  next  cell  and  time  step.  This  procedure  is  the  simplest  of 
“reduced  DSMC”  particle  methods  and  is  known  as  the  “eDSMC”  method  where  the  ‘e” 
stands  for  equilibrium.  (10)  This  approach  was  developed  to  try  to  extend  the  range  of 
stagnation  pressures  for  which  DSMC  could  be  used  to  model  MEMS  and  small  meso- 
thruster  size  flows.  The  approach  was  successful  in  evaluating  the  nozzle  thrust  for  a 
variety  of  MEMS  shapes  and  stagnation  conditions  (11),  however,  it  proved  insufficiently 
accurate  for  modeling  regions  of  the  flow  in  the  vicinity  of  the  nozzle  wall  where  sizeable 
viscosity  exists.  (12)  For  this  reason  we  proceeded  to  examine  the  more  accurate  particle 
techniques  discussed  below. 


2.2  BGK  Approaches 

Although  we  tend  to  develop  computational  techniques  by  artificially  choosing  problems 
where  we  know  they  will  work,  reality  forces  us  back  to  the  situation  where  ultimately  no 
single  approach  will  be  feasible.  In  particular,  micro-  and  meso-  nozzle  transitional  flows 
exhibit  strong  thermo-chemical  non-equilibrium  with  regions  that  are  very  continuum 
and  others  that  are  very  rarefied.  Different  investigators  have  recognized  this  problem 
as  well  and  have  proposed  uncoupled  NS/DSMC  approaches  where  the  NS  equations 
were  used  to  model  the  high  density  portions  in  the  nozzle  and  the  DSMC  was  used  to 
predict  the  rarefied  plume-atmospheric  interaction.  (13)  This  approach  while  demonstrated 
to  be  useful,  has  the  draw  back  that  a  starting  surface  must  be  created  in  the  supersonic 
region  and  the  geometry  should  be  well  defined  in  order  to  keep  the  generation  of  the 
starting  surface  tractable.  For  condensing  flows,  it  requires  one  to  make  the  further 
assumption  that  the  supersaturation  ratio  remains  below  unity  upstream  of  the  starting 
surface.  This  condition  is  not  certainly  not  guaranteed  for  all  stagnation  pressure  and 
nozzle  conditions.  Since  our  goal  is  to  “hybridize”  DSMC  and  BGK  approaches  we  have 
emphasized  the  particle-based  formalism  of  the  BGK  approach  rather  than  finite  volume 
formulations  .(14) 

The  BGK  method  is  a  technique  that  approximates  the  solution  of  the  Boltzmann 
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2.2  BGK  Approaches 


equation, 

8  8  ^8  r°°  r4w 

—  (n/) +  u  •  —  (n/) +  F  •  —  (n/)  =  y  J  n2  (/* ft  -  f  fx)  vradCLdv  (2) 

where  /  and  fi  are  the  pre-collision  single  particle  velocity  distribution  functions  of  the 
two  colliding  particles  and  /*  and  f*  are  their  post-collision  velocity  distribution  functions; 
F  is  an  external  force  per  unit  mass  applied  to  the  particles  (assumed  to  be  zero  for  the 
present  study) ;  a  is  the  differential  cross-section  of  the  binary  collision  and  fl  is  the  solid 
angle;  vr  is  the  relative  velocity  of  the  colliding  particles,  and  n  is  the  number  density. 

To  understand  the  BGK  model,  let  us  concentrate  on  the  collision  term  on  the  right 
hand  side  of  Eq.  2  that  is  difficult  to  compute  because  it  involves  multiple  integrations. 
Different  simplified  models  have  been  proposed  to  model  this  complicated  collision  term 
with  one  such  simplified  model  expressing  the  collision  term  in  a  relaxation  form  known 
as  the  Bhatnagar-Gross-Krook  (BGK)  model: 

4  (nf)  =  "n(fe  -  /)  (3) 

J  collision 

where  n  is  the  number  density,  v  is  the  characteristic  relaxation  frequency  and  fe  is  the 
Maxwellian  distribution  function.  Thus  the  BGK  model  is  a  simplified  linear  approxima¬ 
tion  to  the  nonlinear  collision  term  of  the  Boltzmann  equation  and  provides  an  alternate 
procedure  to  account  for  the  collisional  process  driving  a  flow  toward  equilibrium  without 
modeling  individual  collisions.  This  simplification  is  reasonable  for  flows  with  large  colli¬ 
sion  rates  since  the  details  of  the  particle  interactions  are  not  significant  in  reproducing 
macroscopic  quantities  such  as  temperature,  pressure  or  flow  velocity(15).  Because  the 
BGK  equation  reproduces  correct  moments  and  satisfies  the  H-theorem  for  entropy  pro¬ 
duction,  it  is  accepted  as  an  accurate  physical  model  of  high  density  flows  that  we  will 
discuss  here. 

To  compute  the  collision  term  in  the  BGK  model,  we  first  need  to  define  the  rate 
at  which  the  flow  relaxes  toward  equilibrium.  For  each  species,  the  relaxation  rate  is 
governed  by  the  respective  characteristic  relaxation  frequency  given  by 

v  =  Pr  • nk  Tt~U  (4) 

where  Pr  is  the  Prandtl  number  of  each  species  (assumed  to  be  unity  for  the  BGK  equa¬ 
tion),  k  is  the  Boltzmann  constant,  Tt  is  the  local  translational  temperature  in  a  cell,  and 
/iref  is  the  gas  dynamic  viscosity  of  each  species  at  Tref. 

In  each  cell,  the  number  of  particles  selected  for  each  species  for  the  velocity  sam¬ 
pling  from  the  local  Maxwellian  distribution  depends  on  the  relaxation  frequency  of  the 
corresponding  species  and  the  time  step,  and  given  by 

Nc  =  int(AT(l  —  exp(— vAt)  (5) 

where  At  is  the  time  step,  N  is  the  number  of  particles  of  each  species  in  a  cell  and  the 
int  operator  means  the  nearest  smaller  integer. 
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The  particles,  which  have  been  selected  for  velocity  sampling,  are  assigned  new  thermal 
velocities  according  to  the  Maxwellian  distribution  at  the  local  cell  based  temperature. 


=  —  ln(i?2)  ■ 

■  y/2  kTt/m 

(6) 

=  sin(27ri?i)\/—  ln(i?2)  • 

y/2  kTt/m 

(7) 

=  cos(27ri?3)\/—  ln(i?4)  • 

■  y/2  kTt/m 

(8) 

where  R\  through  P\  are  random  numbers  uniformly  distributed  between  0  and  1  and  the 
superscript  1  indicates  the  velocity  after  the  assignment  of  new  thermal  velocities.  The 
velocities  of  particles  which  have  not  been  preselected  remain  unchanged  in  the  current 
time  step. 

Holway(16)  and  Cercignani(17)  proposed  a  modification  to  the  BGK  equation,  so  as  to 
produce  the  correct  transport  coefficients  of  heat  transfer  and  viscosity.  The  equilibrium 
distribution  function  in  the  ES-BGK  model  is  the  local  anisotropic  three-dimensional 
Gaussian,  fa,  referred  to  as  the  Ellipsoidal  Statistical  (ES)  distribution,  and  the  collision 
term  is  approximated  as 

4  ( nf )  =  un(fG  -  f )  (9) 

J  collision 

The  sampled  velocities  from  the  BGK  procedure  are  modified  to  new  values  conforming 
to  the  ES  distribution  as  follows: 

(10) 


vi  =  Sij  ■  v) 


where  v2  designates  the  modified  velocity  components  and  SLj  is  given  by  the  following 
equation: 


- 


1  —  Pr 
2  Pr 


1 


N 


kTt  N  — 


-  (<  mViVj  >a  <  mVi  >a<  mvj  >a  /  <  m  >a)  -  % 


(11) 

where  the  symbol  <>a  represents  an  averaging  over  all  the  particles  in  a  cell,  StJ  is  the 
Kronecker  delta,  m  is  the  mass  of  a  particle,  and  indices  i,j  refer  to  the  three  Cartesian 
components  of  velocity. 

Similarly,  particles  of  each  species  involved  in  the  flow  field  are  selected  for  the  rota¬ 
tional  and  vibrational  energy  relaxation.  The  corresponding  relaxation  frequency  is  given 
by  the  following  equation 

F, 


1  coll  1  coll  /ia\ 

C  =  uv  =  -z-  (12) 

where  Fcou  is  the  collision  frequency  for  every  species  involved,  and  given  by  the  following 
equation 


Fcoll 


n 


/  kTref 

'  Tt  ' 

V  2tt 

Tref  _ 

'y  1  y  v  ninjAjj  (2  s^) 

i= 1  j= 1 


(13) 


NS[ 


where  n  is  the  overall  number  density,  n  —  ^ 


n7- 


2=1 


7T 


Aij  —  (dref,i  +  drefj) 


I  mnl  +  uij 
rriimj 


(14) 
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In  Eq.  12,  Zr  is  computed  individually  for  every  species  using  the  modified  Parker  (18;  19) 
formula. 


_ 3/5  Zr°° _ 

1  +  (V/72  (T,-/r«5)1/2  +  (7T  +  ir2/4)  (T* /Teq) 


(15) 


where  Teq  is  the  equilibrium  temperature  given  by 


T  = 

-J-  pn 


m 


FrTr  +  FVTV 


eq 


Ft  +  Fr 


(16) 


where  Ft,  Fr  and  Fv  are  the  degrees  of  freedom  for  translational,  rotational  and  vibra¬ 
tional  modes,  and  Tt,  Tr  and  Tv  are  translational,  rotational  and  vibrational  temperatures 
respectively. 

The  constants  T*  and  Zr°°  involved  in  Eq.  15  are  referred  to  as  characteristic  temper¬ 
ature  and  maximum  rotational  collision  number  of  a  polyatomic  species.  For  CO2,  the 
values  of  T*  and  Zr°°  used  in  the  present  work  are  91.5  K  and  18.5  respectively.  The 
vibrational  collision  number  for  CO2  is  assumed  to  be  related  to  the  rotational  collision 
number  by  the  following  relation,  Zv  =  Zr/5. 

Once  we  have  computed  relaxation  frequency  for  rotational  and  vibrational  equilib¬ 
rium,  we  can  find  out  the  number  of  particles  for  each  species  to  be  selected  for  rotational 
and  vibrational  relaxation  from  the  following  equations  respectively,  which  are  in  nature 
similar  to  the  Eq.  5. 


Nr  =  int(AT(l  —  exp(— vrAt)  (17) 

Nv  =  int(A(l  —  exp(— uvAt)  (18) 

In  the  present  work,  for  rotational  and  vibrational  relaxation,  the  particles  are  selected 
from  those  selected  for  velocity  sampling,  i.e.,  among  N„  particles  of  each  species. 

A  comparative  study  of  the  DSMC  and  statistical  BGK  and  ES-BGK  methods  was 
carried  out  in  Ref.  (12)  by  simulating  a  supersonic  expansion  of  argon  to  vacuum  in  a 
nozzle  for  different  pressure  cases  and  comparing  with  the  DSMC  solution  as  the  bench¬ 
mark.  Statistical  BGK  and  ES-BGK  methods  were  found  to  be  in  good  agreement  with 
the  benchmark  DSMC  method  with  ES-BGK  method  showing  a  better  agreement.  It 
was  also  shown  that  the  statistical  methods  of  solving  the  BGK  equation  were  four  times 
more  numerically  efficient  than  the  deterministic  finite  volume  BGK  method  and  twice  as 
efficient  as  the  DSMC  method.  Another  comparative  study  was  carried  out  in  Ref.  (20) 
to  estimate  the  accuracy  of  the  statistical  BGK  and  ES-BGK  methods  for  a  multi-species 
case  involving  involving  the  modeling  of  reinforced  carbon-carbon  (RCC)  thermal  pro¬ 
tection  system  crack  growth  using  kinetic  BGK  method  for  internal  and  external  flows. 
Statistical  BGK  method  was  found  to  be  in  good  agreement  with  the  DSMC  method  for 
the  external  flows.  For  modeling  the  internal  flows  (with  high  local  pressure)  within  the 
vicinity  of  crack,  only  the  BGK  method  was  applied,  because  of  its  cost  effectiveness  for 
the  high  pressure  flows.  Finally,  as  will  be  discussed  in  later  sections,  our  implementation 
of  the  BGK  method  has  allowed  us  to  extend  our  DSMC-based  condensation  modeling  to 
higher  nozzle  pressure  conditions. 
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2.3  Molecular  Dynamics 

In  contrast  to  the  aforementioned  computational  techniques  molecular  dynamics  (MD)  is 
not  a  statistical  method.  Every  particle  is  assumed  to  represent  a  single  true  physical 
particle  and  the  outcome  of  every  collision  is  based  on  the  potential  energy  curve  in 
contrast  to  DSMC  which  assumes  that  the  outcome  is  statistical,  although  based  on  a 
realistic  cross  section.  The  collision  dynamics  are  purely  deterministic  and  obtained  by 
the  integration  of  Hamilton’s  questions.  (21)  MD  is  not  limited  therefore  to  the  binary 
and  dilute  gas  assumption  so  that  it  can  be  applied  to  a  system  of  particles  in  either 
the  solid,  liquid,  or  gaseous  phase.  Because  the  entire  trajectory  is  modeled  the  typical 
applicability  of  MD  is  for  a  system  of  size  on  the  order  of  a  nm  with  106  particles,  and 
for  a  time  scale  of  ~  100  ps,  ,  i.e.,  systems  much  smaller  in  scope  that  are  considered  by 
DSMC  or  its  related  particle  approximations.  In  our  work  we  have  used  MD  as  a  method 
of  modeling  the  “fine-grained”  physics  of  the  problem,  i.e.,  to  obtain  state-specific  reaction 
cross  sections,  improvements  to  the  VHS  total  cross  section,  and  develop  gas-cluster  phase 
collision  and  interaction  models  that  may  be  used  in  “course’-grained  DSMC  and  related 
particle  approaches. 


3  Modeling  of  MEMS  Thrusters 

With  almost  nine  orders  of  magnitude  difference  between  conventional  and  MEMS  nozzle 
thrusters  there  is  no  reason  to  expect  that  the  physics  of  the  latter  should  be  the  same 
as  the  former.  In  fact,  it  is  not.  For  the  MEMS  devices  that  will  be  discussed  in  this 
paper,  the  nozzle  geometries  are  both  rectangular  and  axi-symmetric.  Figure  1(a)  shows 
a  comparison  of  the  stream  wise  velocity,  u/u( 0),  normalized  by  the  value  along  the 
nozzle  axis,  across  the  exit  of  a  3-D  MEMS  nozzle  at  the  nozzle  exit  with  a  conventional 
nozzle.  The  combination  of  the  geometry  and  the  highly  viscous  flow  creates  a  boundary 
layer  which  can  be  seen  to  substantially  envelop  the  core  flow.  The  surface  affects  due 
to  wall  temperature,  heat  transfer,  and  heat  flux  are  substantial  and  will  be  shown  to 
be  controlling  factors  for  micronozzle  performance  and  system  design.  The  right  hand 
portion  of  the  figure  shows  the  development  of  the  boundary  layer  in  a  three-dimensional 
micronozzle  for  a  cold  gas  flow  at  different  axial  locations  from  the  nozzle  throat.  It 
can  be  seen  that  at  the  exit  the  thermal  boundary  layer  engulfs  the  entire  flow  for  these 
conditions.  Yet  fabrication  of  such  systems  in  silicon  and  other  materials  is  well  advanced 
and  concepts  such  as  those  proposed  by  Reed  and  co-authors(22)  shown  in  Fig.  2  have  been 
achieved.  The  figure  shows  an  array  of  thrusters  fabricated  on  a  silicon  chip  with  different 
thruster  sizes  (throat  areas,  nozzle  area  ratios)  available  for  a  range  of  thrust  levels  (from 
/x- N  to  mN)  within  the  same  array.  Several  thrusters  can  be  fired  simultaneously  for  thrust 
levels  higher  than  the  basic  units,  or  in  a  rapid  sequence  in  order  to  provide  gradual  but 
steady  low-g  acceleration.  Decomposing  solid  propellant  pellets  would  be  loaded  into  each 
thruster  prior  to  bonding  the  half  sections.  Decomposing  solids,  also  referred  to  as  gas 
generators,  are  slow-burning  propellants  designed  for  maximum  gas  yield.  Initiation  of 
the  propellants  would  be  accomplished  with  a  diode  laser-based,  fiber-optic  network.  A 
critical  portion  of  our  research  activity  was  the  use  of  DSMC  to  model  the  nozzle  design 
because  it  was  recognized  that  it  will  be  critical  to  the  performance  of  the  system. 
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3.1  Cold  Gas  Thrusters 


Figure  1:  Comparison  of  traditional  and  MEMS  nozzle  boundary  layer  growth. 


The  DSMC  simulations  were  performed  to  understand  in  detail  the  interaction  of 
nozzle-shape,  wall  temperature,  and  stagnation  pressure  and  temperature  effects  on  nozzle 
efficiency,  specific  impulse,  and  thrust.  We  start  with  the  modeling  of  cold  gas  thrusters 
considering  the  affect  of  three  vs.  the  approximate  two-dimensional  model.  Then  we 
consider  whether  better  nozzle  performance  is  obtained  for  higher  temperature  stagna¬ 
tion  conditions.  Finally,  we  examine  high-temperature  gas  thrusters  with  variable  wall 
temperatures  by  use  of  coupled  DSMC/finite-element  method  (FEM)  heat  transfer  sim¬ 
ulations  considering  both  the  affect  of  heat  transfer  on  the  thruster  performance  as  well 
as  the  material  structural  survivability. 

Results  from  various  DSMC  simulations  will  be  presented  below.  Depending  on  the 
stagnation  pressure,  typical  DSMC  simulations  require  about  2.8  to  15  million  particles, 
with  time  steps  on  the  order  of  1  nsec  to  reach  convergence.  Due  to  the  large  computational 
cost  to  perform  a  simulation  of  the  true  three-dimensional  structure,  approximate  two- 
dimensional  and  axisymmetric  structures  were  considered  and  compared  with  the  3-D 
results  to  understand  the  implications  of  using  these  approximate  geometries. 


3.1  Cold  Gas  Thrusters 

Figure  3  shows  the  three  nozzle  geometries  that  were  considered  in  our  studies.  The  left 
hand  portion  of  the  figure  shows  a  three-dimensional  nozzle  of  a  15-deg  expansion  angle 
in  the  XY  plane,  a  throat  width  and  height  of  300  /um,  and  an  area  ratio  of  10.  A  two- 
dimensional  approximation  can  be  constructed  from  this  geometry  by  neglecting  the  end 
surfaces  of  the  nozzle  and  assuming  that  the  height  or  separation  of  the  side  walls  is  very 
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large.  Similarly,  an  axisymmetric  geometry  (RHS  of  Fig.  3)  was  also  considered  keeping 
the  expansion  angle  the  same  as  in  the  two  and  three-dimensional  nozzle  geometries  but 
using  a  throat  radius  of  150  /rm  with  an  exit  to  throat  area  ratio  of  100.  The  test  gas  was 
molecular  nitrogen  with  stagnation  temperature  and  pressures  of  300  K  and  10  kPA,  a  wall 
temperature  of  300  K,  a  Kn=.005  and  a  Re  of  200.  Figure  4  shows  a  summary  of  some  of 
the  important  results  from  our  previous  studies.  (23)  The  comparison  of  the  3-D  versus  the 
approximate  2-D  geometry  for  gas  density  and  temperature  shows  that  the  2-D  geometry 
is  not  a  good  approximation  for  the  real  3-D  configuration.  The  2-D  flow  exhibits  a  higher 
expansion  and  thinner  thermal  boundary  layer  due  to  less  heat  transfer  to  the  nozzle  walls 
compared  to  the  3-D  geometry.  A  similar  conclusion  is  reached  when  comparing  the  3-D 
versus  the  axisymmetric  temperature  and  velocity  profiles.  The  surface  area  to  volume 
ratio  is  less  for  the  axisymmetric  case  than  the  3-D  geometry,  and,  again  the  temperature 
drop  and  velocity  increase  of  the  axisymmetric  case  are  more  consistent  with  isentropic 
supersonic  expansions.  In  contrast,  the  3-D  geometry  does  not  have  a  maximum  velocity 
at  the  exit  and  instead  of  the  temperature  monotonically  falling  it  actually  increases  as 
the  flow  approaches  the  exit.  Further  examination  of  thruster  performance(23)  showed 
that  the  additional  wall  effects  in  the  3-D  geometry  reduce  thrust  by  about  20%  compared 
to  the  2-D  case. 


Figure  3:  Cold  gas  thruster  geometries  (a)  3-D  flat  nozzle  (b)  Axisymmetric  nozzle.  (22) 


3.2  Heated  and  High  Temperature  Gas  Thrusters 

In  this  portion  of  the  research  we  examined  the  sensitivity  of  MEMS  thruster  flows,  with 
the  same  or  similar  geometry  from  the  previous  subsection,  to  increases  in  the  stagnation 
temperature.  As  we  increase  the  stagnation  temperature  we  expect  that  the  thrust  would 
increase  for  the  same  wall  nozzle.  However,  the  increase  in  thrust  is  not  similar  to  that 
which  is  observed  in  traditional  propulsion  systems.  Starting  with  the  velocity  profiles 
for  two  different  stagnation  conditions  of  300  versus  1,000  K,  Fig.  5(a)  shows  that  the 
velocity  profiles  change,  but,  for  both  temperatures  the  profiles  are  dominated  by  surface 
interactions  for  a  3-D  nozzle  with  full  wall  accommodation,  a  throat  Reynolds  number  of 
about  200  and  a  wall  temperature  of  300  K.  The  specific  impulse  for  these  two  cases  is  57 
and  61  s  for  stagnation  temperatures,  Ta,  of  300  and  1,000  K,  respectively.  Although  there 
is  an  increase  in  specific  impulse  it  is  much  less  than  would  be  observed  for  the  comparable 
axisymmetric  cases  because  the  3-D  nozzle  has  a  larger  surface-area  to  volume  ratio.  Any 
factors  which  lead  to  an  increase  in  the  thermal  boundary  layer,  mitigate  the  utility  of 
the  nozzle  to  further  accelerate  the  flow  and  therefore  produce  thrust. 
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3.3  Coupled  Gas-Material  Simulations 


Figure  4:  Cold  gas  thruster  results.  Portions  (a)  and  (b)  show  a  comparison  of  the  flows 
for  a  3-D  (top)  and  a  2-D  (bottom)  geometry  for  gas  density  (kg/m3)  and  temperature 
(K),  respectively.  Portions  (c)  and  (d)  show  a  comparison  of  the  temperature  profiles 
and  stream  wise  velocity,  respectively,  along  the  nozzle  axis  from  the  throat  to  a  region 
downstream  of  the  exit  for  the  3-D  versus  an  axisymmetric  geometry. 


This  can  be  seen  very  clearly  in  Fig.  5(b).  To  reduce  the  computational  effort  (at  that 
time)  the  ratio  of  the  thrust  to  the  thrust  at  the  nozzle  throat  is  shown  as  a  function  of 
distance,  normalized  by  the  throat  radius,  for  an  axisymmetric  nozzle  at  three  different 
stagnation  temperatures,  assuming  a  diffuse,  isothermal  wall  at  300  K.  The  figure  shows 
that  as  the  stagnation  temperature  is  increased  the  nozzle  loses  its  ability  to  accelerate  the 
flow  and  one  should  probably  consider  a  shorter  nozzle  at  higher  stagnation  temperatures 
for  operation  in  the  flow  regime  of  Reynolds  numbers  of  about  200  to  400.  A  similar 
trend  (23)  was  observed  in  other  calculations  that  we  performed  for  the  same  axisymmetric 
nozzle  but  varying  the  stagnation  pressure  for  a  stagnation  temperature  of  2300  K.  We 
found  that  the  lower  pressure  (he.,  more  viscous  flow)  has  a  peak  Isp  closer  to  the  nozzle 
throat,  again  suggesting  that  a  shorter  nozzle  would  be  more  useful. 

3.3  Coupled  Gas-Material  Simulations 

As  the  above  results  and  others  presented  in  Ref.  (23)  show,  the  nozzle  wall  temperature 
and  heat  fluxes  are  major  factors  that  influence  the  gas  dynamics  and  micro-thruster  per¬ 
formance.  The  actual  burn  time  of  a  thruster  is  certainly  an  important  design  parameter 
that  determines  the  impulse  bit  that  will  be  available  for  a  propulsion  maneuver.  In  ad¬ 
dition,  heating  the  microstructure  structure  beyond  its  melting  point  will  certainly  limit 
its  operation.  For  this  reason  we  developed  an  approach  to  simulate  the  coupled  thermal 
and  gas  dynamics.  Since  the  Biot  number  is  sufficiently  small  for  a  microscale  device  the 
lumped-capacitance  method  was  used  to  estimate  a  material  response  time  constant  on 
the  order  of  0.01  s.  This  value  is  five  orders  of  magnitude  slower  than  the  typical  resi- 
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Figure  5:  Warm  gas  thruster  results.  Portion  (a)  shows  the  streamwise  velocity  dis¬ 
tribution  along  the  nozzle  axis  for  two  stagnation  temperature  conditions  for  a  three- 
dimensional  nozzle  and  (b)  shows  a  comparison  of  the  thrust  ratio  for  axisymmetric  noz¬ 
zles  at  three  different  temperatures. 


dence  time  of  the  gas  in  the  micro- nozzle,  i.e.,  10-7  s,  so  that  the  two  calculations  may 
be  loosely  coupled.  (24;  25)  The  coupling  between  the  material  thermal  response  and  the 
DSMC  flow  occurs  when  the  DSMC  heat  fluxes  are  used  as  boundary  conditions  for  the 
heat  conduction  problem.  The  heat  conduction  calculation  then  gives  a  new  wall  temper¬ 
ature  which  is  used  in  the  next  DSMC  simulation.  Figure  6(a)  shows  the  finite  element 
domain  and  its  overlap  with  the  micronozzle  gas  dynamic  domain.  The  3-D  mesh  for  the 
3-D  nozzle  used  about  f600  nodes  and  about  7300  tetrahedral  elements. 

Figure  6(b)  and  (c)  shows  the  gas  and  material  temperature  variation  as  a  function  of 
burn  time  for  a  3-D  nozzle  with  a  stagnation  pressure  of  O.f  atm  and  an  initial  wall  tem¬ 
perature  of  300  K  for  two  different  thermal  boundary  conditions.  The  thermally  insulated 
case  corresponds  to  zero  convective  and  conductive  heat  fluxes.  The  material  (silicon)  and 
the  thruster  boundary  layer  are  isolated.  The  second  condition  of  active  cooling  simulates 
the  presence  of  laminar  water  cooling  in  microchannels  (26)  around  the  perimeter  of  the 
rectangular  material  shape.  Comparison  of  the  two  sequences  of  simulations  shows  that 
in  both  cases  the  nozzle  material  temperature  increases  with  time,  but,  cooling  the  outer 
surface  of  the  thruster  allows  the  material  to  remain  at  a  temperature  that  is  lower  than 
the  melting  temperature.  Therefore  such  a  micro-thruster  would  achieve  a  longer  thruster 
burn  time.  The  gas  flowfields  change  as  well  with  time.  In  the  insulated  case  the  thermal 
boundary  layer  growth  is  higher  than  in  the  cooling  case,  although,  it  is  difficult  from  the 
figures  to  discern  this  trend  quantitatively. 

Instead,  the  difference  in  the  growth  of  the  thermal  boundary  layer  for  the  two  thermal 
boundary  conditions  can  be  better  seen  in  terms  of  the  change  in  thrust  and  coefficient  of 
mass  discharge  shown  in  Fig.  7.  It  can  be  seen  that  for  either  thermal  boundary  condition 
the  thrust  level  falls  as  the  material  temperature  rises,  however,  in  the  cooling  case  the 
material  reaches  a  steady  state  condition  with  the  drop  in  thrust  only  about  6%  compared 
to  the  larger  drop  in  thrust  of  ~  15%  for  the  insulated  case.  Also  shown  is  the  change 
in  the  mass  discharge  coefficient  as  a  function  of  time.  Unlike  the  thrust,  this  coefficient 
is  mainly  affected  by  the  growth  of  the  boundary  layer.  The  more  the  boundary  layer 
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Figure  6:  Coupled  thermal  and  gas  dynamic  computations,  (a)  cross  section  of  compu¬ 
tational  domains  and  time  variation  of  temperature  and  material  responses  for  the  3-D 
nozzle  and  a  stagnation  pressure  of  0.1  atm  for  a  thermally  insulated  boundary  condition 
(b)  and  active  cooling  (c). 


grows,  the  smaller  is  the  effective  cross  sectional  area  of  the  nozzle  and  hence  the  discharge 
coefficient  decreases.  For  an  insulated  nozzle  the  drop  is  as  much  as  55%  of  the  initial 
value  compared  with  only  about  8%  for  the  cooling  case.  The  amount  of  flow  energy  that 
is  lost  to  the  wall  in  the  form  of  heat  can  be  estimated  by  examining  the  total  heat  fluxes 
(sum  of  kinetic  and  internal  energy)  through  a  cross  section  perpendicular  to  the  nozzle 
axis.  For  the  cooling  boundary  condition,  a  non-adiabatic  case,  the  energy  flux  is  not 
conserved  and  decreases  along  the  nozzle  axis.  Figure  4.20  of  Ref.  (9)  shows  that  as  the 
material  heats  up  less  energy  can  be  transferred  to  the  wall  which  in  turn  reduces  the 
decrease  in  thrust  at  later  burn  times.  The  main  reason  that  the  thrust  decreases  in  time 
is  due  to  the  decrease  in  the  mass  discharge  coefficient.  These  results  show  that  the  large 
temporal  variation  in  thrust  and  mass  discharge  coefficient  means  that  the  coupling  of  the 
gas  flow  and  material  response  has  to  be  taken  into  account  in  serious  micro-propulsion 
design  efforts. 


4  Condensation  in  Meso-scale  Expansions 

In  this  section  we  turn  to  the  second  important  physical  aspect  in  small  thruster  expan¬ 
sions,  namely  condensation.  Our  focus  now  turns  away  from  the  physics  of  the  viscous 
boundary  layer  and  is  instead  more  on  understanding  the  extreme  supersonic  expansion 
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Figure  7:  Influence  of  material  response  on  micronozzle  performance  for  a  3-D  nozzle. 
Normalized  thrust  and  coefficient  of  mass  discharge  a  function  of  time. 


and  its  affect  on  condensation.  The  effect  of  the  boundary  layer  is  of  course  still  there  and 
plays  a  critical  role  in  modeling  the  nozzle  lip  region.  Only  a  kinetic  method  can  correctly 
model  the  mass  flux  in  the  back  flow  which  entrains  condensates  that  can  coalescence  and 
stick  to  a  cold  spacecraft  surface. 

There  are  a  number  of  important  challenges  in  modeling  homogeneous  condensation. 
The  first  and  foremost  are  the  multiple  length  scales  in  the  problem  due  to  expansion 
of  the  high  pressure  gas  from  the  plenum  to  the  space  near-vacuum.  Figure  8(a)  shows 
the  three  order  of  magnitude  variation  in  the  Knudsen  number  from  the  nozzle  throat 
to  beyond  the  exit  as  a  function  of  distance  from  the  throat.  As  the  flow  expands  the 
pressure  drops  and  the  supersaturation  ratio,  the  ratio  of  the  saturation  pressure  to  the 
gas  pressure,  increases  by  orders  of  magnitude,  with  the  axial  variation  depending  on 
the  specific  chemical  species,  as  shown  in  Fig.  8(b).  DSMC  cannot  model  the  entire 
flow  and  the  notional  sub-divisions  where  other  approaches  might  be  used  are  shown 
below  as  the  flow  expands  for  an  argon  flow  through  an  orifice.  To  cover  the  wide  range 
of  physical  phenomena  one  needs  both  coarse  grained  and  fine-grained  computational 
techniques.  DSMC  is  a  logical  choice  for  the  coarse  grained  approach,  but,  the  question 
still  remains  about  how  we  can  handle  the  high  pressure  region.  As  will  be  shown  we 
have  successfully  used  the  BGK  approach  to  model  the  plume  from  the  nozzle  exit  to 
the  farfield.  Most  likely  future  hybrid/DSMC  methods  will  open  up  more  opportunities 
to  model  non-laboratory  condensating  flows.  However,  regardless  of  the  choice  of  coarse¬ 
grained  gas  dynamic  approach,  the  appropriate  physical  condensation  models  are  required. 
(Note  that  their  specific  implementation  will  be  different,  however.)  We  have  sought  to  use 
fine-grained  approaches  such  as  molecular  dynamics  to  provide  mechanisms  and  models 
for  cluster  processes  off-line  as  input  for  DSMC  and  BGK  approaches.  However,  the  MD 
results  will  only  be  as  good  as  the  inter-molecular  potential  so  that  the  question  arises 
as  to  whether  there  are  adequate  potential  energy  surfaces  for  even  small,  highly  polar 
polyatomic  systems.  Finally,  abundant  data  exists  to  validate  condensating  flows  from 
small  orifices  and  nozzles,  but,  much  of  it  emphasizes  terminal  cluster  sizes  and  number 
densities  only.  (27;  28) 

In  this  section  we  briefly  outline  the  important  processes  in  condensating  flows  and  the 
two  approaches  that  we  have  used  for  developing  condensation  models,  classical  nucleation 
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Figure  8:  Challenges  of  modeling  multi-scale  condensing  flows,  (a)  Knudsen  number  varies 
by  at  least  three  orders  of  magnitude  causing  the  saturation  ratio  (b)  to  increase  rapidly 
to  greater  than  unity. 


theory  (CNT)  and  MD.  Then  we  provide  selected  examples  of  how  we  have  used  these 
models  in  both  DSMC  and  BGK  simulations.  In  particular  we  will  discuss  our  modeling 
results  of  carbon  dioxide  condensation  and  comparisons  with  the  experiments  of  Ramos 
et  al( 29).  This  comparison  with  this  data  is  important  because  it  provides  a  test  of 
both  the  nozzle  near  and  far-field  cluster  size  and  density  distributions  obtained  from 
the  simulations.  We  then  consider  a  heterogeneous  condensation  case  where  gases  of 
molecular  nitrogen  and  carbon  dioxide  with  different  saturation  ratios  are  both  present  in 
an  expanding  jet  plume.  Molecular  dynamics  was  used  to  calculate  the  sticking  rate  and 
the  simulations  are  compared  with  Rayleigh  scattering  results.  Finally,  we  end  with  an 
example  of  simulating  ethanol  data  where  we  show  that  the  use  of  the  BGK  approach  is 
critical  in  enabling  direct  comparison  with  the  data  of  Wegener  et  al.  (30)  Other  systems 
that  we  have  examined  extensively  include  water(31;  32)  and  ammonia.  (33)  The  ammonia 
system  was  interesting  because  the  condensation  occurs  inside  the  nozzle  and  molecular 
dynamics  was  used  to  check  the  fidelity  of  an  empirical  coalescence  model.  The  cluster  size 
distributions  obtained  from  the  DSMC  simulations  were  found  to  be  in  good  agreement 
with  the  data  of  Bobbert  et  al.  (34) 

4.1  Condensation  Processes 

First  and  foremost  it  should  be  mentioned  that  condensation  changes  the  plume  flowfield 
distribution.  Figure  9  shows  the  DSMC  predicted  plume  temperature  and  mass  fraction 
of  water  gas  for  the  Progress  main  engine  rocket  plume  with  and  without  the  modeling  of 
water  condensation  at  300  km.  (35)  The  Progress  Main  Engine  conditions  were  modeled 
assuming  a  stagnation  pressure  of  approximately  350  kPa,  a  stagnation  temperature  of 
3800  K  and  a  water  species  mole  fraction  at  the  nozzle  exit  of  28.95.%  Figure  9(a)  shows 
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that  the  condensation  process  increases  the  plume  gas  temperature  The  latent  heat  gen¬ 
erated  during  the  condensation  process  is  released  into  the  plume  gases  and  increases  its 
temperature.  Figure  9(b)  shows  that  downstream  of  the  onset  of  condensation,  about 
3.6  m  on  the  centerline,  the  total  gas  mass  density  in  the  condensation  plume  is  less  than 
that  in  the  plume  modeled  without  condensation.  The  reason  for  the  decreased  total  gas 
mass  density  is  that  a  large  fraction  of  water  molecules  condense  into  clusters. 


Mden  (kg/m 
12  1.5e-4 

11  1.1  e-4 

10  9.0e-5 

9  7.0e-5 

8  5.0e-5 

7  3.0e-5 

6  2.0e-5 

5  1.0e-5 

4  3.8e-6 

3  1.2e-6 

2  2.4e-7 

1  7.6e-9 


Figure  9:  Comparison  of  temperature  (a)  and  plume  gas  mass  density  (b)  fields  for  a 
Soyuz  rocket  plume  simulated  with  condensation  (top)  and  without  water  condensation 
(bottom). 


The  processes  that  occur  in  a  condensing  flow  include  nucleation  or  creation  of  clusters, 
evaporation,  sticking  (cluster- monomer  collisions),  and  coalescence  (cluster-cluster)  colli¬ 
sions.  In  the  modeling  of  condensation  processes  in  expanding  supersonic  jets  the  fidelity 
of  the  nucleation  process  is  key.  The  easiest  and  most  pervasive  formalism  to  describe 
cluster  nucleation  is  known  as  classical  nucleation  theory  (CNT).(36;  37)  The  main  idea 
behind  CNT  is  that  clusters  smaller  than  a  critical  size  (based  on  thermodynamic  param¬ 
eters)  will  be  unstable  and  that  clusters  larger  than  the  critical  size  will  spontaneously 
grow  in  a  supersaturated  environment.  The  critical  size,  r*,  is  given  in  CNT  as, 

2o-  c  P 
r*  =  -  ,s  =  — 

pRT In  S'  Ps 


where  a  is  the  surface  tension,  R  is  the  universal  gas  constant,  T  is  the  gas  temperature, 
Ps  is  the  saturation  pressure,  and  S  is  the  saturation  ratio.  Other  important  relationships 
from  CNT  for  the  condensation  (or  sticking),  the  evaporation,  and  the  nucleation  rate  are 
given  respectively  as, 


C 


Tt/I  ,  HJJ 


—  ■  =  — 

yJpinnkTg  ’  \l?2jvmkTc 


pnicr 


(20) 

where  pi  is  the  density  of  the  liquid  state,  ps  is  the  density  of  the  monomer  or  gaseous  state, 
Tg  and  Tc  are  the  temperatures  of  the  gas  and  cluster,  and  q  is  the  sticking  probability.  As 
the  above  relationships  show,  if  one  knows  the  thermodynamic  properties  such  as  surface 
tension,  saturation  pressure,  and  phase  temperatures,  the  above  rates  are  all  defined.  The 
problem  is  that  implicit  in  the  assumptions  of  the  above  rates  is  that  the  volume  of  gas  has 
had  a  sufficient  number  of  collisions  to  equilibrate  (a  long  residence  time) ,  a  condition  that 
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is  certainly  not  expected  in  many  supersonic  expansions  to  vacuum.  Furthermore,  if  the 
conditions  are  such  that  the  critical  cluster  size  is  small,  it  is  not  clear  that  the  macroscopic 
surface  tension  can  actually  describe  a  small  cluster  such  as  a  dimer,  trimer,  etc.  In  fact 
the  CNT  rate  that  was  used  in  the  aforementioned  water  condensation  simulation  of  the 
Soyuz  plume  used  the  correction  to  the  CNT  nucleation  rate  of  Wolk  et  al.  Although  the 
correction  is  motivated  by  experimental  nozzle  measurements  it  is  problematic  because 
for  the  conditions  of  the  Soyuz  plume  simulation  it  varies  by  orders  of  magnitude  and  not 
in  a  consistent  manner. 

Even  considering  a  gas  simpler  than  water  such  as  argon  does  improve  the  uncertainty 
involved  in  using  the  CNT  nucleation  rate.  A  full  kinetic  nucleation  model  for  argon  was 
implemented  and  incorporated  into  the  DSMC.(38;  39)  Figure  10  shows  a  comparison 
of  the  average  cluster  size  (a)  and  the  average  cluster  number  density  (b)  predicted  for 
the  kinetic  and  CNT  models  for  argon  condensation.  The  change  in  the  physics  between 
the  two  nucleation  models  is  evident  from  the  strong  different  in  the  spatial  dependence. 
The  average  cluster  size  is  normalized  to  50  and  1,000  for  the  kinetic  and  CNT  models, 
respectively.  The  average  cluster  number  density  is  normalized  to  8  x  1020  and  8  x  1017  m-3, 
for  the  kinetic  and  CNT  models,  respectively.  The  cluster  number  density  contours  for 
the  simulation  using  the  kinetic  nucleation  process  is  smoother  than  for  the  simulation 
using  CNT  because  computation  temperature  during  a  DSMC  calculation  is  an  inherently 
noisy  process.  In  addition,  the  cluster  number  density  predicted  by  the  kinetic  nucleation 
model  is  about  three  to  four  orders  larger  than  for  CNT. 


Figure  10:  Portions  (a)  and  (b)  show  the  argon  average  cluster  number  density  and  average 
cluster  size  for  two  different  nucleation  models,  respectively. 

This  comparison  shows  that  whenever  possible  it  is  important  to  verify  the  fidelity  of 
CNT  and  when  good  potential  energy  surfaces  exists,  MD  is  a  good  option.  Such  studies 
were  carried  out  by  Li  et  al  to  study  water  condensation  in  small  jets  using  MD  to  obtain 
the  nucleation  mechanisms  and  rates  and  the  sticking  parameters  as  well.  (31;  32)  As 
stagnation  pressure  conditions  increase,  however,  CNT  becomes  more  of  a  viable  option. 

4.2  Carbon  Dioxide  Condensation 

We  present  here  our  modeling  of  carbon  dioxide  condensation  for  the  experimental  con¬ 
ditions  of  of  Ramos  et  al( 29).  In  Ref.  (40)  the  details  of  the  pipe-shaped  nozzle  CFD 
modeling  and  the  specific  starting  surface  for  the  particle  simulations  are  presented.  A 
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comprehensive  DSMC  based  condensation  model  taking  into  account  the  processes  of 
nucleation,  cluster-monomer  sticking  and  non-sticking  collision,  and  cluster  evaporation 
was  developed  and  compared  with  the  experimental  data(29)  in  terms  of  cluster  size  and 
number  density  distributions.  The  experimental  data  was  obtained  by  Ramos  et  al( 29)  as 
part  of  their  quantitative  study  of  cluster  growth  in  free-jet  expansions  of  CO2  for  stagna¬ 
tion  temperature  of  294  K  and  stagnation  pressures  ranging  from  I  to  5  bars.  They  used 
Rayleigh  scattering  and  linear  Raman  spectroscopy  techniques  to  obtain  high  spatial  res¬ 
olution  of  quantitative  growth  of  clusters  in  supersonic  jets  of  CO2.  The  basic  quantities 
measured  along  the  jet  were  the  intensities  of  Rayleigh  scattering  and  of  rotational  and 
vibrational  Raman  lines  from  which  cluster  number  density  and  cluster  size  distributions 
were  derived.  The  work  carried  out  in  Ref.  (40)  involved  simulation  of  condensation  for 
stagnation  pressures  varying  from  1  to  3  bars,  and  the  comparison  with  the  measured 
cluster  number  density  and  average  size  is  shown  in  Fig.  11(a)  and  (b).  It  can  be  seen  that 
the  agreement  between  DSMC  and  experiment  is  reasonable  for  a  stagnation  pressure  of 
1  bar  but  is  unsatisfactory  for  the  higher  stagnation  pressures  of  2  and  3  bar.  DSMC  nu¬ 
merical  parameter  requirements  for  higher  stagnation  pressures,  such  as  4  and  5  bar,  could 
not  be  satisfied.  Further  comparison  of  the  DSMC  simulated  carbon  dioxide  translational 
temperature  (which  was  found  to  be  equilibrated  with  the  rotational  temperature)  pre¬ 
sented  in  the  left  portion  of  Fig.  12  shows  that  the  near-field  agreement  corresponding  to 
the  onset  of  condensation  is  not  well  described  by  the  simulations  for  stagnation  pressures 
of  1  and  3  bar. 

The  inability  to  model  the  entire  stagnation  pressure  range  of  data  and  the  suspicion 
that  the  modeling  of  the  release  of  heat  due  to  condensation  might  not  be  correct  led 
to  our  recent  efforts  to  model  condensation  directly  in  BGK.(41)  In  the  first  phase,  the 
statistical  BGK  scheme  was  implemented  in  the  baseline  DSMC  code  for  multi-species 
and  polyatomic  gas  flows  using  Eqs.  9  through  18.  In  the  second  phase,  a  new  heat  ad¬ 
dition/accommodation  model  was  incorporated  in  the  BGK  and  DSMC  schemes.  The 
general  concept  behind  the  heat  addition/accommodation  model  is  that  the  added  heat 
scales  up  the  random  velocities  of  particles  by  a  factor  which  is  determined  by  the  ratio 
of  heat  added  per  cell  per  time  step  and  the  total  random  translational  energy  of  all  the 
particles  belonging  to  the  cell.  The  local  translational  temperature  of  the  cell  is  then 
computed  based  on  the  scaled  up  particle  velocities.  The  rotational  and  vibrational  ener¬ 
gies  also  become  scaled  up  because  molecular  collisions  transfer  translational  to  rotational 
and  vibrational  energies.  The  remaining  portions  of  Figs.  II  ((c)  and  (d))  and  12  (center 
and  right  hand  side)  show  the  improvement  of  the  new  heat  addition  model.  The  cluster 
number  density,  average  cluster  size,  and  the  gas  temperature  close  to  the  condensation 
onset  near  the  nozzle  exit  are  now  in  good  agreement  with  the  experimental  data.  The 
results  showing  the  new  heat  accommodation  model  use  the  BGK  simulation,  but,  the 
same  results  (not  shown  here)  for  the  one  and  three  bar  stagnation  pressure  cases  were 
obtained  with  the  DSMC  approach  as  well.  What  was  not  possible  to  simulate  with  the 
DSMC  model  is  the  5  bar  case  and  as  Fig.  12  and  Ref.  (41)  show  the  simulations  nicely 
capture  the  trends  for  all  measured  profiles  as  a  function  of  stagnation  pressure. 

Having  compared  our  DSMC  and  BGK  simulations  for  carbon  dioxide  condensation 
with  experiment  we  wanted  to  see  whether  we  would  be  able  to  model  a  mixture  of  a 
carbon  dioxide  and  molecular  nitrogen  expanding  flow.  Figure  8(b)  shows  that  the  su- 
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Figure  11:  Comparison  of  cluster  number  density  (a)  and  (c)  and  average  cluster  size  (b) 
and  (d)  between  simulations  and  experiment.  The  top  row  shows  the  results  of  DSMC 
simulations  using  an  older  heat  addition  model.  The  bottom  row  shows  a  comparison 
of  the  old  versus  the  new  heat  addition  model  and  the  experiment  for  the  three-bar 
stagnation  case  using  BGK. 


Figure  12:  Comparison  of  gas  temperatures.  Left  figure  shows  older  heat  addition  results 
compared  to  experiment,  center  figure  shows  a  comparison  of  older  and  present  results 
using  the  new  heat  addition  model  for  the  3  bar  stagnation  condition,  and  right  figure 
shows  the  comparison  of  data  and  new  heat  addition  model  for  three  stagnation  pressures. 


persaturation  ratio  reaches  unity  for  carbon  dioxide  sooner  than  it  does  for  molecular 
nitrogen,  as  expected  from  common  experience.  Figure  13(a)  shows  the  potential  energy 
curves  for  dimers  of  carbon  dioxide  and  nitrogen.  The  shallower  well  for  nitrogen  com¬ 
pared  to  carbon  dioxide  is  consistent  with  its  later  condensation.  The  saturation  ratio 
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profiles,  S,  profiles  shown  in  the  figure  for  the  two  species  corresponds  to  an  expansion 
from  a  15  degree  conical  nozzle,  with  throat  diameter  of  0.6  in,  stagnation  pressure  of 
1.02  atm,  stagnation  temperature  of  285  K,  and  a  mixture  of  95  N2  and  only  5%  C02.(28) 
Because  the  carbon  dioxide  condenses  prior  to  nitrogen  even  a  small  amount  of  it  will  form 
clusters  to  then  enable  nitrogen  condensation  to  occur.  This  is  borne  out  by  the  Rayleigh 
scattering  measurements  presented  in  Fig.  13(b)  which  shows  that  the  pure  nitrogen  flow 
is  indistinguishable  from  our  DSMC  simulated  noncondensing  flow  but  when  even  a  trace 
amount  of  carbon  dioxide  gas  is  added  condensation  can  now  be  observed. 

We  proposed  that  the  initial  carbon  dioxide  dimers  are  created  by  triple  collisions, 

C02  +  C02  (C02)* 

(C02);  +  N2  (C02)2  +  N2 

where  the  triple  collision  rate  depends  on  the  lifetime  and  collision  frequency  of  the  col¬ 
lision  complex  (C02)$.  The  dimer  formation  probabilities  for  this  process  were  obtained 
from  MD  simulations  for  temperature  conditions  that  correspond  to  the  nucleation  re¬ 
gion  in  the  plume.  Since  the  stagnation  pressure  was  sufficiently  low,  DSMC  simulations 
were  performed  using  generalized  cluster-monomer  non-sticking  collision,  condensation, 
and  coalescence  models. (40)  Figure  13(b)  shows  that  the  DSMC  simulation  that  included 
this  condensation  mechanism  (pink  line)  compares  well  with  the  Rayleigh  scattering  mea¬ 
surements  (green  dots). 


.......  Noncondensation  Flow 

A  95%  Nj  +  5%  C02,  DSMC  C02Homo. 

♦  100%  N2,  AEDC  Data 

•  95%  N2  +  5%  C02,  AEDC  Data 
..........  95%  n2  +  5%  C02,  DSMC  Heter. 
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Figure  13:  Comparison  of  Lennard-Jones  potential  energy  curves  for  carbon  dioxide  and 
molecular  nitrogen  (a)  and  comparison  of  Rayleigh  scattering  measurements  with  DSMC 
simulations  of  heterogeneous  carbon  dioxide-nitrogen  condensating  flows  (b). 


4.3  Ethanol  Condensation 

Ethanol  jet  expansions  have  been  studied  by  Yarygin  et  al  as  a  surrogate  for  hydrazine 
back  flow  contamination.  (42)  In  these  experiments  they  used  model  thrusters  in  a  vacuum 
chamber  to  test  the  effectiveness  of  screens  in  blocking  the  contaminants  from  reaching 
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the  backflow  area  of  the  plume.  Such  experiments,  while  relevant  to  operational  con¬ 
ditions,  are  difficult  to  simulate.  Instead,  we  turned  to  the  measurements  of  Wegener 
et  al.  (30)  who  carried  out  experiments  to  study  homogeneous  nucleation  of  ethanol  for 
varying  amounts  of  ethanol  concentration  into  air  at  the  same  total  mixture  pressure. 
The  conditions  for  the  cases  that  we  modeled  are  given  in  Table  1.  In  our  first  attempt 
at  simulating  these  experiments(43)  we  found  that  it  was  not  possible  to  use  DSMC  to 
simulate  the  actual  experimental  conditions.  Therefore  an  approximation  was  made  in 
the  simulations  by  assuming  that  the  total  mixture  pressure  was  that  of  the  pure  ethanol 
pressure  actually  used  in  the  experiments,  i.e.,  in  Ref.  (43)  it  was  assumed  that  the  total 
pressure  was  1.87  kPa  instead  of  the  true  total  pressure  of  83.46  kPa.  In  our  more  recent 
work(44)  we  were  able  to  model  the  actual  experimental  conditions  using  the  statisti¬ 
cal  BGK  approach  with  a  condensation  models  reformulated  to  take  advantage  of  the 
BGK  approach.  The  condensation  models  included  CNT-based  nucleation  and  evapora¬ 
tion  models  as  well  as  coalescence  and  sticking.  In  addition  a  new  weighting  scheme  was 
implemented  to  reduce  the  computational  load. 

Figure  14  shows  the  important  results  of  the  simulations.  Part  (a)  shows  a  comparison 
of  our  simulated  cluster  to  ethanol  gas  mass  fraction  with  the  experiment  for  four  different 
ethanol  mass  fractions  as  a  function  of  distance  from  the  nozzle  exit.  Compared  to  the 
previous  work(43)  (seen  in  (b))  the  improvement  is  quite  good.  The  supersaturation 
pressure  for  ethanol  as  a  function  of  temperature  is  shown  in  part  (c)  to  provide  a  reference 
for  the  experimental  and  simulated  condensation  onset  value.  In  part  (d)  we  compare  our 
prediction  on  the  condensation  onset  using  the  original  definition  of  Wegener  et  al  as  the 
location  in  the  flow  where  the  condensed  mass  fraction  reaches  10-4.  Both  simulated  and 
measured  condensation  onset  occur  in  the  flow  at  ethanol  gas  pressures  and  temperatures 
above  the  saturation  pressure,  as  expected.  Although  the  data  of  Wegener  et  al  is  old, 
it  is  consistent  with  the  more  recent  data  of  Wyslouzil  et  al( 45).  Our  numerical  results 
show  good  agreement  with  the  experiments,  thus  confirming  the  use  of  the  BGK  based 
condensation  model  for  computationally  challenging  high  pressure  flow  applications. 


Table  1:  Wegener  et  al( 30)  Experimental  Conditions® 


Case 

Ethanol  Supply 
Mass  Fraction 

Ethanol  Partial 

Pressure  (Pa) 

1 

0.008 

211.4 

2 

0.005 

163.5 

3 

0.0034 

89.7 

4 

0.0025 

65.7 

“Stagnation  pressure  and  temperature  are 
83.4  KPa  and  296  K,  respectively  for  each  case. 
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Figure  14:  Comparison  of  simulations  and  measurements  of  the  cluster  to  ethanol  mass 
fraction  for  different  supply  mass  fractions  of  ethanol  (a)  and  our  previous  modeling  (b). 
The  saturation  pressure  for  ethanol  is  shown  in  (c)  and  a  comparison  of  our  predicted 
onset  point  with  data  and  previous  work  is  presented  in  (d). 

5  Future  Directions 

The  improvements  in  modeling  MEMS  and  meso-sized  propulsion  systems  in  terms  of 
thrust  as  well  as  contamination  has  seen  dramatic  improvements  over  the  last  decade.  At 
the  same  time,  MEMS  fabrication  techniques  have  continued  to  improve  leading  to  higher 
stagnation  pressures  while  maintaining  small  sizes.  The  higher  stagnation  pressures  have 
the  desirable  effect  of  increasing  thrust,  as  expected,  and  in  addition  in  reducing  the  large 
viscous  boundary  layers  at  lower  stagnation  values.  While  much  has  been  fabricated  and 
demonstrated  in  a  laboratory  environment,  the  non- academic  environment  has  been  slow 
to  embrace  micro-propulsion  even  for  small  satellites.  Perhaps  the  only  systems  actively 
being  considered  today  are  colloidal  micro-thrusters. 

Contamination  concerns  on  large  satellites  are  generally  handled  operationally  by  siz¬ 
ing  and  building  in  extra  capacity  to  meet  the  mission  lifetime  requirements.  Of  course 
such  a  strategy  will  not  work  for  resource  limited  constellations  of  small  satellites.  More¬ 
over,  as  other  forms  of  propulsion  such  as  electric  and  mixtures  of  electric  and  chemical  are 
considered  it  can  be  expected  that  the  spacecraft  contamination  environment  will  be  even 
more  challenging.  New  propellants  consisting  of  ionic  liquids  condensing  on  spacecraft  sur¬ 
faces  in  the  atomic  oxygen  low  earth  orbit  environment  will  present  plenty  of  challenges 
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to  future  missions.  In  terms  of  modeling  and  simulation,  we  are  only  starting  to  obtain  in 
inkling  of  the  new  chemical  models  that  will  be  necessary  to  model  contaminating  sprays 
of  these  advanced  propellants. 
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